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Chapter-! 


1 A General Introduction with Historical Perspective 

1.1 Role of Computational Fluid Dynamics in Engineering Research 

Computational Fluid Dynamics has attained a vital role in 
fluid mechanics research to the level of having its own standing 
based on analytical and experimental knowledge of the flow physics 
as well as the tools of numerical analysis. Analytical solutions, 
albeit being derived for simplified forms of governing equations, 
have the distinct advantage of immediately highlighting some of 
the fundamental features of the flow physics. However, real life 
situations are usually so complicated that analytical solutions 
are not available in most of the cases and experiments are too 
expensive. Due to the rapid growth of efficient and economic 
computing facilities as well as numerical modelling methods, the 
predictive procedures support experiments, enrich and even extend 
the range of analytical solutions and finally contribute actively 
to product development. 

1.1.1 Applications of Computational Fluid Dynamics 

Major applications of numerical simulation are related to 
aerospace research, combustion studies, nuclear reactor control, 
mechanics of reacting fluid flow, meteorological predictions etc. 
In these applications, the rapid decrease in the cost of 



computations ( yfith the advent of high speed modern-day computers) 
compared to the staggering cost of performing state of the art 
experiments, has made numerical simulation an economical 
alternative. Efficient numerical algorithms are being used and 
continuously upgraded to substitute rigorous experiments on wind 
tunnel testing or combustion research. However, the numerically 
predicted models may further need fine-tuning through final 
experimental tests. Besides economic advantage, numerical 
simulation has the potential to depict the detailed flow field 
information accurately. It can even monitor flow conditions 
involving fast transients, high temperatures, toxic substances 
etc., which pose considerable difficulties in experimental 
measurements . 

1.1.2 Salient Features of Numerical Simulation Techniques 

Modelling of real life phenomena starts with the basic 
conservation principles of physical quantities like mass, 
momentum, energy and species. These fundamental principles are 
often expressed in terms of algebraic, ordinary or partial 
differential, or even integral forms of equations. These governing 
equations along with appropriate boundary and/or initial 
conditions are mathematically reformatted into equivalent discrete 
representations of the problem by suitable spatial and/or temporal 
discretization techniques. The discretized equations are, in turn, 
solved by appropriate solution procedures. Finally, the numerical 
solutions thus obtained, are processed further to calculate 



quantities of interest 

It may be mentioned that the present wide-spread acceptance of 
the computational methods of simulation as an essential dimension 
to the studies on mechanics of fluids is the outcome of effective 
numerical predictions in a variety of applications. In all such 
studies accuracy of the numerical data largely depends on the 
correctness of the mathematical approximations in the modelling of 
the physical situation. This is especially true for turbulent flow 
predictions, combustion, modelling of reacting fluids, nuclear 
reactor studies etc. for which rigorous mathematical descriptions 
are not yet available. However, even for problems in which the 
mathematical model is well understood, dealing with irregular and 
complex geometries and selection of the final numerical solution 
scheme are as much important as the mathematical modelling. In all 
these cases, very many important and interesting phenomena are 
often revealed by numerical modelling. 


1,2 Complexities Associated with Numerical Modelling of Fluid 
Flow Problems : 


Incorporation of proper flow physics is a very crucial step in 
the numerical prediction process since the simulation by itself 
cannot reproduce physics that has not been duly considered. 
Physical complexities that are often encountered in many flow 
situations are compressibility, chaos and turbulence, fast 


are 



transients in the flow field etc These are usually taken care by 
suitable constitutive relations or modifications in the basic 
Navier Stokes and mass conservation equations. Often, heat and 
mass transfer also occur simultaneously in the flow situations. 
Many of these problems finally lead to a set of coupled partial 
differential equations which need to be solved simultaneously. 
However it is possible to decouple variables like enthalpy or 
temperature from the solution of the basic velocity and pressure 
fields in some situations. 

While solution of a complicated real life problem can be 
handled by a set of add-on relations with the basic model of 
Navier Stokes equations, geometry of the flow domain necessarily 
requires careful treatment. An orthogonal coordinate frame and a 
matching set of orthogonal velocity components are always a 
perfect choice which provide clarity in the adaptation of any 
numerical scheme, in the solution procedure and during the final 
flow visualization. However, in real life problems, alignment of 
the flow variable layouts with respect to a fixed standard 
coordinate frame may not prevail. Also, the velocity/pressure 
boundary conditions may be available only at a set of points lying 
on an irregular boundary so that adoption of a standard coordinate 
frame turns out to be complicated and often impossible. Many 
approaches have been tried in the past to handle complex 
geometries, such as (i) rectangular, axisymmetric or polar 
coordinates with step wise approximation, of curved boundaries ( so 


called castellated grids), (ii)u3age of orthogonal curvilinear 
coordinate systems, (iii) development of general nonorthogonal 
body conforming coordinates etc. Each of these approaches has 
specific usefulness to the respective numerical scheme. 

The selection of a numerical scheme for estimating the flow 
physics is a critical step. This is because of the fact that even 
the incompressible constant property viscous flows involve 
nonlinear partial differential equations. The situations are much 
more complicated in presence of any of the above mentioned 
physical or geometrical complexities. Aspects which require closer 
scrutiny in a numerical scheme for incompressible viscous flows 
are the proper accounting of the nonlinear convective kinematics 
and the pressure field evaluation. The convection kinematics needs 
incorporation of the upstream biased modelling, especially, for 
convection dominated flow situations. Also, in unsteady flows, 
excepting the explicit schemes, most of the others need iterative 
velocity evaluation because of nonlinear convection terms. 
However, getting an acceptable pressure field is a difficult task 
since there is no direct equation for evaluating pressure, in 
incompressible viscous flow. The non-zero divergence of an 
incorrect velocity field is often employed to adjust pressure 
field till cell wise incompressibility of the fluid is achieved in 
all cells. Different schemes have so far been developed to 
optimize this process of preserving incompressibility and 
simultaneously extracting pressure, corresponding to various flow 



s ixua uions 


1 . 3 Aim of the Present Work 

In this thesis, the task of developing a general calculation 
procedure for viscous incompressible flows in complex geometries 
has been undertaken. It has been decided to solve the unsteady 
Navier Stokes eq.uations for both orthogonal and non-orthogonal 
grids and follow a time marching procedure to reach a steady state 
(if the flow field is physically steady) or to reach a particular 
time level (if the flow field is unsteady periodic). The 

importance of continuity equation as a constraint on the velocity 

field has been deployed in solving the pressure field. Finally, 

the aim of the present thesis is to employ this algorithm for 
analyzing the wake-zone-aerodynamics behind partially enclosed 
bluff bodies. 

1.4 Development of Computational Methods in Fluid Mechanics 
- A historical Perspective 

1.4.1 Flow Formulations : 

The solution of different forms of Navier Stokes equations has 
been the subject of research for the past three decades. 

Concurrent attempts have been successfully made with different 
sets of variables for describing these momentum balance equations. 



Two dimensional analysis based on the application of stream 
function and vorticity variables, has the advantages of removing 
pressure from the set of solution variables and simultaneously 
satisfying continuity requirements. However, extension of such a 
stream funct ion-vorticity approach to three dimensional flows is 
not possible. Primitive variable formulation involving velocity 
components and pressure, on the other hand, is used by 
comparatively a large number of researchers since this is directly 
extendable to three dimensional geometries. Such a formulation 
often meets the requirements of prescribed boundary conditions of 
application problems, in quite a straight forward manner.lt may be 
mentioned that in finite element and other related approaches, the 
interpolation functions of the stream function and vorticity must 
have at least continuous first derivatives needing Cl class 
element. The primitive variable formulation, however, needs only 
CO class element, although satisfying zero divergence of 
velocities over each element sometimes poses serious challenges. 
1.4.2 Methods on Regular Geometries 

As mentioned earlier, solution of Navier Stokes equations with 
primitive variables has two basic challenges namely, (i) absence 
of any obvious equation for solution of pressure and, (ii) the 
non-linearity of the convection kinematics. Also, the nature of 
the equation is such that the velocity and pressure solutions may 
exhibit non-physical oscillations (wiggles) if a collocated grid 
is used for finite difference formulation of all the variables. 


However, this well known problem of ’red -- black checker board 
splitting’ of pressure is circumvented through either the 
implementation of a staggered grid arrangement or with the use of 
special interpolation techniques on collocated regular grids. 

(i) Primitive Variable Based Solvers: 

The full Navier Stokes equations in conservative form have been 
solved by finite difference based explicit transient algorithms 
such as MAC and SMAC by Harlow and Welch[l] and Harlow and 
Amsden[2] respectively. In these algorithms the provisional 
velocities at any instant of time, are first estimated from 
momentum equations using pressure and velocity values of the 
earlier time step. These estimated velocities, however, are yet to 
satisfy continuity requirements. If the cell wise divergence 
values of these estimated velocity fields are non-zero, then there 
must be some amount of mass accumulation or annihilation in each 
cell which is not physically possible. This indicates that the 
pressure field is incorrect. The non-zero divergence of velocity 
can be utilized to correct pressure and the velocities in turn, 
and a Poisson type pressure correction equation can be obtained. 
The procedures of MAC and SMAC differ in the implementations of 
this pressure velocity iteration to arrive at zero divergence. 
These methods use a layer of imaginary cells around the boundary 
of the physical domain necessitating update of boundary conditions 
after every change in internal velocity and pressure values. In a 
related technique of artificial compressibility method developed 


by Chorin[3], simultaneous iteration on pressure and velocity 
components are performed to obtain a converged solution 
Viecelli[4l has shown the equivalence of all the above mentioned 
methods. Hirt and Cook[5] have developed a simplified explicit 
flow solver by making use of the aforesaid pressure-velocity 
iteration for a wide variety of problems involving three 
dimensional incompressible flows. While these explicit schemes are 
quite efficient in the studies of temporal flow development, they 
have stability restrictions on the time step values. Implicit 
methods are indeed attractive since they have no such restrictions 
regarding the choice of time steps. The finite volume method based 
SIMPLE algorithm of Patankar and Spalding[6] provides an efficient 
implicit procedure. Enhancement of this basic SIMPLE code has been 
made by Patankar himself through the version of SIMPLER[7,8] and 
by Vandoormaal and Raithhy through SIMPLEC[ 9 ] . In SIMPLER, at every 
iteration level or time step, a first estimate of velocity is 
obtained from the evaluation of the convection diffusion parts of 
the momentum equations. The remaining source terms and the 
pressure terms constitute a pressure equation. The pressure and 
velocity are iteratively corrected in much the same way as in 
SIMPLE. Although the changes necessary to incorporate SIMPLEC into 
a SIMPLE code are minor, the consequences can be great as SIMPLEC 
eliminates the approximations made in SIMPLE while deriving the 
pressure velocity corrections. Jang, Jetli and AcharyaClO] 
reported a comparative illustration between the operator splitting 


PISO[ll] algorithm of Issa and the SIMPLE family Of algorithms 
(ii) Convection Diffusion Modelling : 

Numerical modelling of the convection diffusion phenomena 
demands as much consideration as deriving a correct pressure 
solution. Most of the above algorithms and works are quite stable 
in this regard and have good convergence characteristics 
corresponding to their ranges of applications. Just as algorithms 
were generated aiming at proper solution of pressure, parallel 
developments occurred for the efficient handling of nonlinear 
convection terms in the flow equations. It is well known that 
numerical instabilities arise when the flow equations contain 
dominant first derivative terms. Initial convection discretisation 
models were space-centred second order difference formulae. 
Instabilities with such formulations were characterized by 
oscillations (wiggles) in the velocity solutions. It was found 
that these wiggles could be eliminated by using one-sided (upwind) 
first order finite difference approximations for the derivatives 
in the convection terms. Solutions for intermediate Reynolds 
number flows using such upwind approximations have been published 
by several investigators[ 14 , 15 , 23 ] . However, based on detailed 
studies on the nature of these differencing schemes, Roache[ 16 , 17] 
and several others [ 18 , 1 9 ] have critisized the idea of introducing 
’artificial viscosity’ into the physical problem through upwind 
schemes. DeAllen et al [20], Spalding[21] Raithby et al[22] and 
others [23] have modified the upstream differencing concept to 



blend its advantages with space centred differencing Patankar[7], 
in his hybrid scheme, has made use of different approximations 
including this first order differencing to closely follow the one 
dimensional analytical solution of steady state equation. Many 
successful usages of similar multiple approximations based on 
cell Reynolds numbers have been reported in literature 
[10-14,24-27]. Han et al [39] have used a hybrid differencing 
scheme in order to solve the flow problem in a lid— driven square 

5 

cavity for a Reynolds number of 1x10 . In the context of high 
Reynolds number turbulent flows, the upwinding discretization 
schemes due to Leschziner and Rodi [40] are particularly notable. 
Roache[16] proved that the truncation error generated due to the 
use of first order accurate difference formula may become too 
large to be ignored. He gave a mathematical estimate of this error 
and showed that as compared to the second order accuracy of all 
other terms, this error may be fixed with the corresponding 
truncated term of the Taylor series interpolation of convection. 
The coefficient of the second derivative in the truncated 
convection interpolation formula has been identified to introduce 
a numerical diffusivity for the transported quantity, A similar 
truncation error arises out of backward first order time 
differencing in transient studies. As the steady state is 
approached, this error approaches zero. Restricting the Courant 
number less than unity, however, takes care of this problem. 
Raithby[28] has given a host of critical examples to emphasize the 



applicability and limitation of first order npwinding Depending 
on the nature of the problem, several upstream-weighted 
formulations have been attempted in between. Skew upstream 
dif f erenc ing [ 29 ] of Raithby, skewed positive influence coefficient 
upwind procedure of Schneider and Raw[30] , mass-flow-weighted 
skew upwind procedure due to Hassan et al [31] and the controlled 
numerical diffusion with internal feedback (CONDIF) scheme of 
Runchal [32] are a few important efforts in modelling of 
convection terms. 

Another new concept which has been established by Galpin and 
Raithbyt33] is the treatment of the non-linearities in the 
incompressible Navier Stokes equations by Newton-Raphson 
linearization . Higher order upstream differencing of convection 
terms had also been studied widely by Pollard and Siu[34], Davis 
and Mallinson[ 35 ] , Leonard [36-38], Han et al [39], Leschziner and 
Rodi [40], Shyy [41] and Raithby et al [42]. The second order 
upwind schemes, especially QUICK of Leonard[38], and SUDS of 
Raithby [29] have received a lot of appreciation for their 
effective applicability and stability. Tao and Sparrow[43] and the 
series of communications from Patel et al [44-45] have discussed 
in detail about the performance characteristics of all the so far 
available upwind schemes. Vanka[46] has presented a comparison 
among the second order upwind schemes of Leonard {QUICK}[49], 
Wilkes and Thompson[47] and Atias et al[48]. Gresho and Lee have 
demonstrated that the optimal upwind difference scheme with an 



effective grid Peclet number of hyperbolic tan function, is 
unconditionally stable Barrett [50] in his Super upwind scheme has 
given an estimate of distribution of transported perturbation on 
both sides of a point of original disturbance. 

Zeinkiewica and co-workers have successfully used upstream 
asymmetric basis functions [ 51 ] and modified function based 
upwind[52] in corresponding Finite Element formulations. Huyakorn, 
Taylor, Lee and Gresho[53] compared different mixed interpolation 
Finite Element methods for fluid flow problems. Sani et al[54] 
have conducted an exhaustive study on spurious pressure field 
generated in certain FEM analyses of fluid flow problems and also 
outlined other aspects of associated difficulties. A higher order 
convection model in FEM has been suggested by Donea and 
Giulliani [ 55 ] . Kato et al[56] reported to have solved three 
dimensional cavity flow for large Reynolds numbers using GSMAC FEM 
algorithm which they have considered optimum on the basis of 
computational load. 

1.4.3 Methods on Complex Geometries : 

Finite difference and finite volume based methods have been 
widely applied to problems in regular geometry whose boundaries 
fall on the coordinate lines of a simple orthogonal coordinate 
system. The orthogonality of the coordinate frame makes it 
possible to decouple the integrands so as to carry out individual 
line integration products for evaluating the surface integrals. 
Thus, many researchers have tried to generate orthogonal grids for 



complex solution domains while more complex geometries are 
necessarily modelled with general non-orthogonal curvilinear 
grids , 

(i) Orthogonal Grids: 

Usually in engineering problems such regular orthogonal grids 
may be unsuitable or an analytic coordinate frame may not be 
available. In order to extend the aforesaid established flow 
solvers to the complex geometries, use of generalised orthogonal 
body-conforming coordinate systems have been made in different 
combinations. While the mathematical utility of such orthogonal 
grids in fluid mechanical problems is very attarctive, 
availability of a smooth grid is equally restrictive. Hung and 
Brown[57] developed an implicit MAC like two-dimensional scheme 
using grid-oriented orthogonal set of velocity components. 
However, their formulation was based on non-conservative form of 
equations. Pope [58] applied the conservative forms of the 
equations with finite volume method to compute turbulent 
recirculating flows in diffusers. Gosman and Rapley[59] used a 
similar procedure for fully developed flows in ducts of arbitrary 
cross-section. Turbulent flow calculations were done on such grids 
by Rapley[60], Habib and Whitelaw[61] and Rastogi[62]. Lawal and 
Majumdar[63] have performed calculations in the entrance region of 
an arbitrary shaped cross section. Assuming a parabolic flow, 
collocated cartesian velocity components on orthogonal grids in 
each cross-sectional plane have been solved along with a Poisson 


equation for pressure 
(ii) Non Orthogonal Grid 

In non-orthogonal systems, general curvilinear coordinates are 
employed which exactly coincide with the boundaries of the 
physical domain. The interior also has smooth grids in terms of 
slope and curvature continuity. The governing partial differential 
equations are transformed in terms of these curvilinear 
coordinates and finite difference representations are derived in 
the transformed domain, Gal-chen and Somerville [ 64 ] and Faghri et 
al[65] used algebraic coordinate transformation. While Gal-chen 
used a MAC-like grid arrangement coinciding with the grid lines, 
Faghri et al had used covariant components of velocities. Thompson 
and coworkers [66,67] have made remarkable contributions to the 
development of numerical grid generation techniques for solving 
elliptic partial differential equations related to both external 
and internal flow problems. Steger and Sorenson [68,69] developed 
a computer code to solve Poisson equations with grid control near 
the boundary and with coarse-fine sequential solver. Vanka et 
al[23], Maliska and Raithby[70], and Shyy et al[71] have extended 
the available solvers of different flow problems to irregular 
geometries. The basic flow solvers have been based on the family 
of SIMPLE algorithms. The adaptive grid procedure of Acharya and 
Patankar [72] for parabolic flows may be mentioned here for its 
detailed discussion on the choice of different discretization 
schemes and their performance in the adaptive grid layout. A new 



differencing scheme, the w-tracing scheme, had been suggested here 
which showed extremely good matching with the corresponding 
analytical curve. Patel and Briggs[73] used MAC in boundary fitted 
curvilinear coordinates. Karki and Patankar [74] used a 
generalised coordinate system with co-variant velocities for flows 
in complex geometry. In 1981, the first successful application of 
finite volume methods using non-orthogonal coordinates and 
collocated grids were reported by Hsu [75]. Concurrent works had 
also been reported by Prakash[76] and Rhie[77]. Disproving the 
impracticability of collocated arrangements, several other works 
had been regularly reported since then. A few of them are due to 
Rhie and Chow[78], Rhie et al[79], Peric[803, Burns et al[81] and 
Reggio and Camarero [ 82 ] . It may be worthwhile to note from Bernard 
and Thompson [83] that while staggered grid is a direct solution 
to avoid pressure split with conventional finite difference 
interpolations of flow variables, proper interpolation of 
velocities still remains crucial to avoid oscillatory velocity 
field even on a staggered grid. An appropriate pressure 
interpolation avoids pressure splits in these methods of 
collocated grids. Use of collocated velocity grids has the 
following obvious advantages ( Peric et al[84] ): 

(a) all velocity variables share the same location implying 
only one set of velocity control volume; 

(b) the convection contribution to the coefficient of the 
discretized equations is same for all the velocity variables; 



(c) for complex geometries cartesian velocity components can 
be used, in conjunction with non-orthogonal coordinates yielding 
simpler equations than when coordinate oriented velocity 
components are employed; 

(d) there are fewer constraints on the numerical grid, since 
there is no need to evaluate the so called "curvature terms" , and 

(e) treatment of boundary conditions and implementation of 
higher order differencing are simpler* 

Further detailed discussion on these aspects are available in 
the studies of Hsu[75], Peric[80] and Miller and Schmidt[85]. A 
finite volume based procedure has been very successfully developed 
by Majumdar [86] using such collocated velocities and pressure. 
The concept of "momentum interpolation" of cell face pressures 
from nodal values has been established showing equivalence in the 
methods of staggered and collocated arrangements! 87 3 . However, 
such momentum interpolation needs under relaxation between 
iterative levels of velocity pressure adjustments. Majumdar [88] 
has also clarified that a consistent under relaxation will enhance 
the rate of convergence and at the same time achieve a unique 
solution that is independent of the under relaxation parameter 
used. Williams! 89 ] has applied a collocated primitive variable 
formulation on two dimensional cavity problem where he introduced 
an amount of intermediate compressibility into a projection scheme 
to obtain Helmholtz type expression for a pressure variable. 

Solution of the Navier Stokes equations on irregular geometries 



being the state oi-the-art problem, Baliga and Patankar[90j opened 
a new dimension of the SIMPLER algorithm with the successful 
proposition of the control volume based finite element 
method ( CVFEM ) . By discretizing the domain into six noded macro 
elements each of which consisted of four, three-noded triangular 
subelements, the discretized equations were obtained by deriving 
algebraic approximations to integral conservation equations 
applied to these polygonal control volumes. While pressure 
variables were stored only at the vertices of the macro elements, 
all others were - stored at nodes of the subelements. Pressure 
interpolation was linear while other dependent variables were made 
sensitive to cell Peclet number and the direction of the velocity 
vectors, Prakash and Patankar[25, 26] used equal order 
interpolation for pressure and velocity for these CVFEM 
algorithms. Schneider and Raw[30] introduced special upstream 
biased convection diffusion model for CVFEM while Muir and Baliga 
[91] extended them to three dimensional applications. Considerable 
overall improvement has been reported by Prakash[24] by 
introducing special shape functions even for the source terms in 
the transport equations. 

For flows in complex geometries, the finite element method 
appears to be an attractive alternative owing to its intrinsic 
geometrically additive properties. Application of this method to 
fluid dynamical problems had been started with Stokes flow. As 
mentioned earlier, preference is given to primitive variable 



formulation over the vector potential (stream function) formulation 
because of algorithmic suitability. Howeverj as far as concepts of 
variational calculus are adaptable to FEM, the rigorous theory of 
extremizat ion subject to (differential) constraints is available 
only for linear elliptic boundary value problems concerning 
creeping flow situations, Oden[ 92] , Olson[93], Baker[94] and Taylor 
and Hood[95] extended the finite element methods successfully to a 
variety of fluid flow problems. While in structural mechanics 
problems, the mass conservation is intrinsic because of the 
Lagrangian description, for flow problems it needs as much 
attention in FEM as in FDM. Zienkiewicz [96] proposed penalty 
function method for linear systems, viz., Stokes problem . It was 
later extended to two dimensional incompressible flows by 
Reddy[97]. Because of stability and convergence requirements, 
reduced integrations and corresponding variable layouts had been 
prescribed by Zienkiewicz et al[98] and Reddy[99]. Tuann and Olson 
[100] preferred solving the incompressible flow equations in 
primitive variables using continuity equation to generate the 
Poisson equation for pressure. However, modelling of convection 
kinematics remained unresolved until Christie et al [51] used 
asymmetric linear and quadratic basis functions to overcome the 
difficulty of significant first derivatives in a one dimensional 
two-point boundary problem. Heinrich and Zienkiewicz[ 52 ] extended 
this idea to two dimensional applications. A mixed interpolation 
coupled finite element solution algorithm of the velocity and 
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pressure for incompressible flow bss been discussed in detail by 
Taylor and Hughes [101]. Comini and Del Guidice [102] have 
calculated provisional velocity values based on an estimated 
pressure distribution and momentum interpolation by making use of 
a Finite Element technique. They have obtained appropriate 
corrections to the provisional velocities through satisfying 
continuity requirements . Snith and Brebbia [103] described a method 
for the solution of transient, incompressible viscous flow in two 
dimensions. The dependent variables of stream function and 
vorticity, were approximated over each triangular element using 
linear interpolation functions. As an application, a study of the 
vortex street development behind a rectangular obstruction has 
been described. 

In order to use velocity-pressure formulations in FEM, 
appropriate interpolation functions for velocity and pressure 
variables are needed. For error consistency reasons, the inter- 
polating functions for pressure should be polynomials of lower 
degree than those for velocities. Usually, quadratic 
interpolating shape functions for velocities and linear shape 
functions for pressure are deployed. To date, two types of 
mixed-interpolation finite elements have been commonly used. 
These are eight noded serendipity and six-noded triangular 
elements. The first type of elements was recommended by Hood and 
Taylor [104] and the second type of elements are preferred by 
Kawahara et al , [105], However, despite valuable contributions 



from Baker t 106 , 107 ] , Zienkiewicz and Godbole[ 108 ] , and Lohmeyer 
et al [109] finite element algorithms still have some unresolved 
difficulties such as lack of transportive property in the 
advective terms. 

1.5 Scope and Objective of the Present Work : 

Despite several stimulating developments as stated in the 
earlier sections, a satisfactory answer with regard to the 
question of computing incompressible flows in complex geometries 
is still not available. The common difficulties are, namely, 
derivation of the conservation principles for non-orthogonal 
control volumes, development of appropriate interpolation schemes 
and the enormity of book keeping involved with respect to both 
dependent and independent variables. 

In the Present work, a novel algorithm (EXTRA-FLAG) for 
non-orthogonal geometries using primitive variables is proposed 
which resolves some of the above mentioned difficulties. Integral 
forms of conservation equations for finite control volumes in the 
domain are considered. These forms are particularly appealing 
because of the potential for using powerful numerical quadratures. 
Ideas such as element wise interpolation and transformation of 
non-orthogonal element geometry into a square master element which 
are commonly used in FEM, have been incorporated. Collocated and 
Coordinate oriented velocities are used. Pressure nodes are 
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staggered with respect to the velocity nodes 

The proposed EXTRA-FLAG algorithm has been applied to several 
test problems to establish accuracy, stability and applicability. 
The problems are basically unsteady isothermal two dimensional 
incompressible flows in orthogonal and non-orthogonal geometries. 
The name of the algorithm has been chosen as EXTRA-FLAG which 
signifies Exp licit TRansient Algorithm for FLows in Arbitrary 
Geometry. The algorithm has been finally employed to analyze the 
structure of wakes behind partially enclosed bluff bodies. 



23 


Chapter-2 

2. Description of the EXTRA-FLAG Algorithm ; 

This algorithm is based on integral mass and momentum balance 
applied over non-rectangular control volumes in general. The 
various steps involved in the algorithm are as follows. The 
integral flow equations are first converted into algebraic form 
through numerical quadrature, coupled with spatial interpolation 
of the velocity and pressure fields in terms of their nodal 
values. In order to facilitate numerical quadrature, each 
non-orthogonal (quadrilateral) control volume is mapped into a 
standard rectangular cell with the help of local curvilinear 
coordinates. Explicit time-marching is employed for the momentum 
equations using pressure field corresponding to the previous time 
level. The satisfaction of the continuity principle on the other 
hand, is enforced in an implicit fashion by iteratively correcting 
the velocity and pressure fields. Thus, the proposed algorithm 
combines the powerful features of interpolation and numerical 
quadrature used in Finite Element Method with the physical insight 
of the control-volume approach employed in MAC, SIMPLE and related 
algorithms. Due to the adoption of non-rectangular control 
volumes, the present method can be used for solving flow problems 
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in arbitrary-shaped geometries. The solution procedure (called 
EXTRA-FLAG because it is an EXplicit TRansient Algorithm for Flows 
in Arbitrary Geometry) is explained in detail in the following 
sections 

2.1 Derivation of Governing Equations for a Control Volume : 

For the illustration of the proposed algorithm, let us consider 
two-dimensional unsteady incompressible laminar flow inside an 
arbitrary shaped geometry. As shown in Fig.2.1, the domain is 
discretized into small curvilinear quadrilateral cells. The 
velocity nodes are located at the vertices and the pressure nodes 
are located at the centroids of these cells, In order to determine 
the values of the velocity and pressure variables , the momentum 
and continuity principles are applied over the corresponding 
control volumes which surround the velocity and the pressure nodes 
respectively. For instance , the continuity control volumes (CVl) 
are formed by grid lines connecting velocity nodes (Fig. 2. 2a) 
Similarly momentum control volumes (CV2) are formed by curvilinear 
quadrilaterals whose vertices are pressure nodes (Fig. 2. 2b) . It 

is worthwhile to note that in the present algorithm, the velocity 
components are collocated at a velocity node so that the same 
momentum control volume can be employed for both x- and y— 
momentum equations. 

Considering a typical control volume for mass continuity, CVl , 
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the integral form of the mass balance equation can be written as 

a 


^ JJ pdA + Jp V . n = 0 


( 2 , 1 ) 


C V 1 


Csl 


where dZ is an elemental length on the boundary (CSl) of the 
control volume and n is the local outward unit normal vector. For 
incompressible flow, the first term of equation (2-1) is 
identically zero. Therefore, the equation can be simplified as 

( 2 . 2 ) 


r 


t v , nd't = ( u n^+ V n^ ) d^ ~ 0 

1 csi 

where n and n are the direction cosines of the outward normal 

^ y 

n on the boundary CSI. 

Now considering a typical momentum control volume CV2 , the 
integral form of x- and y- momentum balance equations can be 
written as 


3t 


u(pdA) + cpuCpV 


cvi 


and 


€■5 2 


n 


d-C) = F. 


(2.3a) 


at 


v(p dA) + <Pv(pV . nd-C) =F 


cv 2 


CS2 


(2.3b) 


In the above equations , F^and F^ are the components of the 
resultant forces acting on the control volume in x- and y- 
directions respectively. These can be evaluated as follows : 


= i .9 ^ • n d-C 


CS2 


P g dA 


cva 
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p dA (2.6) 

For incompressible flow, the body force contribution in 
equations (2.5) and (2.6) can be absorbed into the corresponding 
pressure term. The final expressions for the x and y momentum 
equations in integral form are ; 

udA + <>u(un + vn)d'& 

X y 

CS2 
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For evaluation of the area and line integrals occurring in the 
above expressions , it is necessary to define interpolation scheme 
for the variation of each flow variable over the concerned control 
volume or control surface. Also, suitable numerical quadratures 
are to be employed to perform integration over curvilinear areas 
and lines. These aspects are described in detail in the ensuing 
sections . 

2.2 Interpolation of Flow Variables Within Respective Control 
Volumes : 

Since the order of the momentum equations is higher than that 
of the continuity equation , the velocity variables have been 
interpolated using bi-quadratic functions while applying the 
momentum balance principle ; on the other hand , for the 
satisfaction of continuity , bilinear interpolation functions have 
been used to describe the velocity variables over the respective 
control volume. Mathematically these interpolation schemes can be 


written as follows. 
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For the continuity equation : 

“k = Ml 


u = 


k«l kVl 

and for momentum equations (2.7) and {2.8) : 


(2.9) 


u = 




V = 




( 2 . 10 ) 


k a 1 k * 1 

The summation in equation (2.9) is performed by making use of 
the velocities at four vertices of the continuity cell (Fig. 2, 2a). 
Similarly the summation in equation (2.10) is accomplished on all 
the neighboring velocity nodes surrounding the momentum control 


volume boundary , as shown in Fig. 2. 3. Also are 


respectively bi-linear and bi-quadratic interpolation functions. 
The expression for these functions are presented in the next 
section. 

The variation of pressure over the momentum control volume is 
also expressed in the form : 


p = 



( 2 . 11 ) 


It should be noted , however , that the above-discussed 
interpolation schemes for the flow variables are based on 
intuitive reasoning ; alternative schemes are definitely possible 
which may provide superior handling of specific terms in the flow 
equations for particular applications. The proposed algorithm is 



Fig. 2.3 Interpolation domain for momentum bale 
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general in the sense that any such interpolation scheme can 
suitably incorporated. 

As regards the transient terms in the momentum equations ^ 
lumping concept has been employed for the sake of simplicity 
this is easily achieved by setting 
r r 


at 


cvi 
r r 


at 


u dA = p A u 

^ CV2 


V dA = p 


( 2.12 ) 


be 

the 

and 


CV2 

where A^^^ is the area of the control volume CV2 and u , v are 
the time derivatives of the velocity components at the central 
node of the control volume , CV2 . 

2*3 Numerical Evaluation of the Area and Line Integrals : 

The integrals appearing in equations (2.7) and (2.8) have been 
evaluated using Gauss-Legendre quadrature , since they are very 
complex to evaluate analytically . The first step required for the 
application of such numerical quadratures is the mapping of the 
non-orthogonal control volume into a master square cell whose 
coordinates vary from (-1) to (+1) in both the directions (Fig. 
2,4). The mapping of the control volume itself can be achieved 
through what is known as the isoparametric transformation. In this 
transformation , the coordinates x and y for any point within the 
non-orthogonal control volume are interpolated in terms of their 
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nodal values in much the same way as the flow variables . When 
quadratic interpolation functions are used , quadratic curves in 
the physical domain will be mapped as straight lines in the 
transformed domain. Therefore a curvilinear control volume bound 
by quadratic curves will be exactly mapped onto a square 
computational cell as shown in Fig. 2.4. The interpolation scheme 
for the coordinates of a momentum control volume (CV2) are given 
by : 


X = 


kVl 


k k 


y = 


kTi 


k k 


(2.13) 


where ( ) are the coordinates of the nodes surrounding 
the control volume as shown in Fig, 2.3. For transforming the 
continuity cells (CVl), linear interpolation is used in the form : 


=1 


M, X, 
k k 


- 1 . ^ 


(2.14) 

_ It ■ K 

k a 1 k » 1 

Equation (2.14) implies that the continuity control volumes are 

quadrilaterals made up of straight boundaries. Such linear 

interpolation of coordinates is adequate since velocity components 

are also interpolated by linear functions for the satisfaction of 

the continuity principle. The coordinate transformation is 

isoparametric in this case as well . 

The standard square cell in the above-mentioned transformations 
is represented in terms of local coordinates 5 and IJ . The range 
of variation of these coordinates over the cell boundary is from 
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(-1) to (+1) . Thus 5 after isoparametric transformation , the 

limits of all the integrals can be conveniently set as (”1) to 
(+1). The interpolation functions and are also neatly 

expressible in terms of the local coordinates S and T] . These are 
given by : 

N^= ~ ? ri (?-2)(r)-2) ; ^ 2 " ^ ^ ( 4-^^ ) (Ti-2 ) ; 

N3= ^ ? T (e+2)(i7-2) ; N^=: ^ 5 (?+2)(4V); 

Ng= ? r] (?+ 2 }(r]+ 2 ) ; ^5= 32 ( 4 - 5 ^ ) (t1+2 ) ; 

N^= “ 5 n (5-2)(T7+2) ; Ng= ^ 5 {?-2)(4-rj^}; (2.15 ) 

^ 9 = (4-5^)(4-r|^) and 

= f- (l-5)(l-r7) ; = -7 {l+?){l-r)) ; 

M 3 = ~ (1+5){1+T1) ; = -7 (l-5)(l+n) . 

For the completion of the transformation process , it is 


necessary to convert the arc lengths , areas , surface integrals , 
derivatives and integrals in the physical coordinate system in 
terms of the computational coordinates 5 and 17 . Such conversions 
for the momentum control volume (CV2) will now be illustrated . 
Applying the chain rule , the relationships between lengths can be 



or 
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written as : 



But , for the faces 1 and 3 , 77 = constant and for the faces 2 and 
4 , 5 = constant (see Fig. 2. 4) . Also , on sides 3 and 4 , the 
elemental length vector is in the direction opposite to the 
increasing ^ or 77 respectively , when an anti-clockwise sense of 
integration is adopted for all the cyclic boundary integrals. In 
order to be consistent with the anticlockwise sense of 
integration, the outward normals at the cell face are considered 
as positive. 


Using the above conditions in equation (2.19) , it can be shown 

that 
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where t , n (s = 1,2, 3, 4 ) are respectively the local unit 
s s 

tangential vector and the unit outward normal vector on the sth 
side . Also , the denominators of the above expressions can be 
calculated as 



The elemental length , d-t , along each side can be represented 


as : 


d-C = 
s 

r at • 
1 a? . 

) d? 

s 

for 

s = 1 and 3 ; 

dl = 
s 

f dt ' 
1 at] , 

) dn , 

s 

for 

s = 2 and 4 ; 


Lastly , the elemental area , dA, appearing in the transient 
terms of the equation (2.7) and (2.8) can be expressed as : 

dA - I J 1 de dn (2.22 ) 

where the | J 1 is the determinant of the Jacobian matrix [ J ] 


40 


defined in equation (2.17) . Now all the integrals appearing in 

equations (2.7) and (2.8) can be rewritten in terms of the 
transformed coordinates 5 and rj . The Gauss-Legendre quadrature 
can be easily adopted for evaluating these integrals since the 
limits are (-1) to (+1). The quadrature formula for the area 


integrals is given by 

+,i +,i 


(2.23 ) 


A - 


F(5,T7) d? drj = 2 ^ 

k » 1 1*1 

where Gauss sampling points and w^ , w^ are the 

corresponding weights . The coordinates of the sampling points 
( known as Gauss points ) are the zeros of the m th degree 

Legendre polynomials , where m is the total number of Gauss points 
used in a particular direction. The line integrals can be 
similarly evaluated as : 


+,i 


F{?) de 


- 1 


k = 1 




(5,) 


W, 


and 


F(ri) dr] 


= 1 = 

1-1 


(n,) 


K, (2.24 ) 


In the present work , 3- point quadrature has been employed for 
evaluating the momentum equations . The corresponding Gauss point 
coordinates and weights are provided in Table 2.1 . 



Gauss Point Location and Weights for 3-point < 


Gauss 

Co-ordinate 

Weight 

point 

location (§ or n) 

1 

- / 0.6 

0 . 5555555556 

2 

0.0 

0.8888888889 

3 

- / 0.6 

0.5555555556 


the application of the above quadrature formu] 
u , V and p , their x- and y- derivative 
cosines n , n of the local outward 

y y 

b of Jacobian | J [ , elemental length df- 

area dA are evaluated at each Gauss poini 
boundary or control volume and all such Gau; 
ons are summed up according to equations (I 
Eventually , matrix equations are obtained 
ilume which represent the x- and y- momentui 
. The solution of the equations is describe 
on . For the continuity equation (2.2) , i 

to employ the numerical quadrature since the 
simpler and the velocity interpolation is lin 
tion of the continuity equation is also discuss 
section. 

X Equations for Momentum and Mass Balance : 
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The final discretized momentum equations in matrix form are 


u 

^ V j 
2x1 


[C] + [D] 


(u) 


+ [ S ] { p } 


. ( V) J CV2 

18x1 2x4 4x1 


C V2 


(2.25 ) 


2x18 

The convective coefficient matrix in the momentum equations can 
be expressed as : 


where 


[ C ] = 

2x18 


(Clj) = (C2j) = 4) 
1x9 1x9 


(Clj) 


(C2j)J 


(2.26a ) 


1x9 1x9 


N f u n + V n ) dl (2.26b ) 

P j X y 


CV2 


Applying Gauss-Legendre quadrature , 


(Cij) = (C2j) 


= — t 


CV 2 s»l k*l 


N. ) (u n + V n ) N 

J At, q X q y q 

q- 1 


w, A-C- 

k S 

(2.26c ) 


where A 


cv 2 


[ D ] = 

2x18 


\ J 1 

w. w 

1 

JC L 

can be 

expressed 


i 1 

(Du;-;r 

r (Dv^/) _ 

1x9 

each 


(2.27a ) 


where 
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Applying Gauss-Legendre quadrature , 
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The pressure matrix also can be expressed in a similar way : 


[ S ] 


■ L J 


2x4 1x4 each 


where 


'"ti* = - WT 


{ s ) = - i 

^ 2j^ PA, 


M n d-t 
J X 


M n d-t 
J y 


(2.28a ) 


{2.28b ) 


Applying Gauss-Legendre quadrature , 




^ " PA 


— [ ! t (M. n ) 

CV2 L S=1 fc=l ^ ^ 

— [ t t (M n ) 

CV2 L s^i k4-i J y k,; 


w, 

k S 


w, 

k s 


{2.28c ) 


The mass balance equation (2) can be rewritten as : 


{un + vn )d'&=0 

X y 


y{un +vn )A^=0 
Zj, s xs s ys s 

s» 1 


(2.29) 


where u^ , are respectively the average values of x and y 

velocity components on the sth side and n , n are the 

X 5 y s 

direction cosines of the outward normal for that side . The 
appropriate evaluation of u^ and v^ on each face of the continuity 
cell is somewhat intricate and this is discussed in the subsequent 


section . 
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2.5 Derivation of Pressure Correction Equation 

As the name of the algorithm (EXTRA-FLAG) suggests , the 
velocities are evaluated explicitly at each momentum control 
volume for the next time step by making use of the momentum 
equations (2.25) . The pressure , however , is evaluated 
implicitly in conjunction with the continuity equation . Equation 
(2.25) can be rewritten as 



where the superscripts n and ntl denote time levels . 


In the implementation of the algorithm , only those velocity 
and pressure values which satisfy both the momentum and the 
continuity equations are given the superscript n+1. The 
provisional values which have been obtained through the explicit 
calculations of equation (2.25) ( and not satisfying continuity ) 
are represented with the superscript * . Thus, equation (2.25) can 
be rewritten as; 



Note that in equation (2.31) , the pressure is still at the nth 
level since it has not been so far updated via the compliance of 
continuity equation. Equation (2.31) can be expressed in a 
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slightly more general form as 


u 


(u) 




.(v) . 


-I- At 


[C] [D] 


(u) 


^(v) 


n 


+ At [ S ] 


{ P*} (2. 


32) 


It is conjectured that after satisfying the continuity equation 
in each cell , provisional values ( values with * superscripts ) 
for the dependent variables will be transformed into final 
converged values for that time step . Thus the accomplishment of 

is performed through an 


u 


n+1 * 

-4 u V ” 


n+ 1 j * 
V and p - 


n+1 


iterative procedure which will be discussed herein , 

Subtracting equation (2.32) from equation (2.30) the expressions 
for the velocity corrections can be obtained as 


n+1 * 

u - u 

n+1 * 

V - V 


= 

ISvj 


S At [ S ] < 5p 


(2.33) 


2 X \ 2x1 2x44x1 

where the 5u , 3v are the velocity corrections and 6p is the 
pressure correction. Equation (2.33) can be used to calculate the 
velocity corrections, if the pressure corrections are known . 

The evaluation of the pressure corrections is carried out 
through an implicit satisfaction of the continuity equation as 
described below. Assuming that the mass balance principle is 
satisfied by the converged values of velocity components at both 
the time levels of n and n+1 , one can write 


s= 1 


n 

( u n 


n 


+ v n ) A^ = 
s xs s ys s 


0 


( 2.34a ) 
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_n+l n+1 

) ( n + V n 

s^i s xs s 


) = 0 
ys ' s 


The unconverged velocities 


however 


( 2.34b ) 
do not satisfy 


(un +vn )A-C.=R 
s xs s ys s 


continuity and therefore leave a non-zero residue . Thus , 

( 2.34c ) 

Subtracting equation (2.34c) from equation (2.34b) , we obtain : 

( 2.35 ) 


I 


I 


( 5u n + 5v n ) = R 

s xs s ys s 


where 6u and dv are the average velocity corrections on the 

o S 

Sth side , so as to satisfy continuity principle in the control 
volume , From equation (2.34c) and {2.34a) » the residue can be 

written as : 


R = 



S“ 1 


( ) n + ( - V ) n 

s s xs s s ys 




= t ( «u 

S» 1 

where 5u 


n + 6v n 1 

mom,s xs mom,s ys j s 


and 6v 


(2.36 ) 

are the average velocity 


mom,s mom,s 

modifications in the sth side of the continuity cell arising from 

the explicit time increment calculation of the momentum equation 

(2.31) . Combining equations (2.35) and (2.36) , 


(6u n +5v n )A-t = / f 

s xs s ys s A.. V 
3=1 s = 1 


du 


mom 


n + 5v n 1 A-f 

, s xs mom,s ys j s 


(2.37 ) 

By substituting for the velocity corrections in the LHS of 
equation (2.37) in terms of pressure corrections , equation (2.37) 
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can. be converted into a pressure correction equation . The crucial 

step in this process , however , is the estimation of the average 

velocity corrections 6u and 5v in terms of the nodal velocities 

s s 

of the sth side . Referring to Fig. 2. 5 , the obvious choice for 

the average which may arise in one’s mind is : 


5u^ = (5uj + 5ujj) / 2 ; 5v ^ ~ + ^v^^) / 2 

/ 2 : = (5vj, + / 2 ) 


> (2.38 ) 


and so on . The nodal velocity corrections Su^ , Sv^ etc. , can in 
turn , be written in terms of the pressure corrections via 
equation (2.33) . It was found that such a simple averaging 

procedure leads to a pressure correction equation which exhibits 
decoupling between the pressures of odd and even continuity cells. 
This pressure splitting has also been observed by Majumdar, Rodi 
and Vanka [87] and it has to be avoided since it introduces 
numerical oscillations in the pressure field . 

The pressure splitting can be avoided by a proper definition of 
the average corrections 5u and 5v on each side . Defining 

3 S 

face-centre velocity corrections 6u„ and 6v„ , the average 

X ) S X I s 

corrections can be calculated as follows • On the east side of CVl 
( refer Fig. 2.5 ) 


5u = 

e 6 1 1 6 


6v = 
e 


X 

6 -ii 


6v „ + — =• 5v 


. -1 

5u „ 

III 3 

f , e 

. -4 

^v „ 

I II 3 

f ,e 


( 2.39 ) 



2.5 Pressure-velocity coupling for a continuity cor 
volume. 
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similarly , the average corrections on the sides ( n , w , s ) can 
also be estimated in terms of the corresponding nodal and 
face-centre values of velocity corrections . Note that equation 
(2.39) amounts to the application of Simpson’s rule along the 
boundaries of the continuity control volume CVl . 

It now remains to link this face-centre velocity corrections 
with the pressure corrections . A close look at equation (2.33) 
provides the key to resolve this is.sue . Equation (2,33) can be 
viewed as an expression relating the correction in the velocity 
vector , 5V to the gradient of pressure correction ''^(Sp). Indeed , 
this simplified momentum equation can be rewritten as : 

5V = - ^ V (6p) (2.40) 


Now , the face-centre velocity corrections on the east face 

r - 5p^- 5p^ 


and 5v„ are given by : 

t f e 


5u„ = ( i . 5v) . 

f , e ” f , e 


or 


5u 


f , e 


P 


At 


E P 


■£ -^P 


5Pp- 


IT 


M 


(2.41a) 


similarly , ^ } (2-41b) 

In the above equations e is the unit vector connecting the 

central node P within the control volume to the pressure node E of 

the eastern neighbouring continuity cell ( see Fig. 2.5 ). Also, 

A^ is the distance between the nodes P and E . In a similar 
E 
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fashion , the face-centre velocity corrections for the north , 
west and south sides of the control volume can be obtained as : 


At ( ^(w,w,s) ~ j 

^f , {n,w, s) ~ p ^ 

(N,W,S) 

At ( ~ ] 

f,(n,w,s)“p 


6p - 6p 

Al 

(N , W, S) 

5p - 6p 

•tD -Er/irti ri\ 


where the subscripts within the brackets may be chosen 


(2.42a) 


(2.42b) 


appropriately for the side under consideration . 

These relationships can be used in the derivation of the 
pressure correction equation in the following way . The mass 


residue contribution R to the east face is given by 

e 

A-C ( 5u n + 5v n } = - R 
e e xe e ye e 

substituting for 5u , 5v from equation (2.39) , 

e e 

AZ -f —4 5u_ + —4 5u_, + du. _ 1 n + 


(2.43a) 


f ,e 

r n 

J xe 

5v^ 

1 " 

f 

,e / 

the 

nodal 


^ 5v + -4- 5v + -| 5v_ In = - R (2.43b) 
e \ 6 II 6 III 3 f,ej ye e 

Now , from equation (2.33) , the nodal velocity corrections are 

given by : 


- At [ S ] 


(2.44a) 


! !, T ■'! 


52 


and 


r 5u 

f III 

I 5v 

L III 


At t S ] 


III 


5p 


NE 


(2,44b) 


J 

as can be seen from Fig. 2.5. In equations {2.44a) and (2.44b) the 
subscripts II and III refer to quantities corresponding to nodes 
II and III . Thus, substituting equations (2.41a), (2.41b}, 
(2.44a), and (2.44b) into equation (2.43b), the residue 
contribution of the east face R can be calculated in terms of 


e 

pressure corrections. Similarly for the other faces ( n , w , s ) 
also, the velocity corrections at nodes and face-centres can be 
written in terms of appropriate pressure corrections. The final 
expression for the pressure correction equation can be then 
derived from the mass residue equation (37), in a form : 


CP 




n + 
xs 


5v 



R (2.45) 

mom 


1X9 9X1 

where CP is the coefficient matrix corresponding to the pressure 
correction array and R is the mass residue generated after the 

mom 

explicit evaluation of the momentum equations through equation 
(2.31). The left hand side of equation (2.45) has a strong 
diagonally dominant structure , with the 6p of the concerned 
continuity cell linked to the pressure corrections of all the 
eight neighbouring cells ( see Fig. 2. 5 ) . 

The pressure correction equations written using equation 
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(2.45) for all the continuity cells can be solved together by 
matrix inversion technique , line-by-line sweeping or by 
point-by-point method . In the present work , point-by-point 
iteration has been employed for the sake of simplicity . It 

may be mentioned that after the explicit momentum evaluations , 
only changes which have taken place in the nodal velocity 
values are available . The changes in the face-centre velocities 
are not available . To start with , these may be estimated by 
'simple averaging as given in equation (2*38) and the approximate 
value of the residue R can be evaluated . As pressure values get 
modified during the continuity iterations , the corresponding 
changes in the nodal or face-centre velocities are also calculated 
using equation (2.33) and (2.41) respectively . The change in the 
mass residues due to the changes in the velocities are again 
evaluated and used in equation (2.45) . In this fashion , the 
pressure and velocity corrections are iteratively improved until 
the corrections to the mass residues decrease below a pre-defined 

* S “7 

small value («# 10 to 10 ) for all the continuity cells . The 

(n+l)th level velocity and pressure values at each node are then 

(2.46 ) 


calculated as : 


n+1 

P 

n 

= P 

+ 

5p 

n+1 

u 

n 

= u 

+ 

5u 

n+1 

V 

n 

= V 

+ 

5v 
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where 5p , 6u and Sv are the overall pressure and velocity 
corrections accumulated iteratively over the time step At * Thus 
by explicit momentum evaluation and implicit satisfaction of 
continuity i time marching of nodal velocities and pressures are 
continued until steady state is reached from given initial 
conditions . 

An important point to be noted is that in the present 
algorithm , the domain boundaries always contain velocity nodes . 
Since known velocity conditions are more often prescribed ( 
through no-slip or given flow velocity ) at the domain boundaries 
, the velocity corrections for the corresponding nodes and 
face-centres can be set equal to zero . Therefore , the 
contributions to mass residue (or correction in the mass residue ) 
will be zero for these boundaries . 


2.6 Non dimensional Form of Governing Equations and Definition of 
Problems Under Investigation 

The equation of continuity (2.2) and the momentum equations 
(2.7) and (2.8) are expressed in nondimensional forms as 

r 

<)(Un+Vn )d'&=0 
X y 


3t 




U dA + o U (U n + V n ) d-t 

y 


CS2 


(2.47) 
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cb ^ - P 


ciz 


X * Re dX "^Re I ay ^ dx J ^yj 


(2.48) 


rr 


dt 


cU 


V dA + J V (U + V riy) d-t 

CS2 


} = 


ii 

n X 2 

av 1 ( 

■ au 

av 

1 \ 

fi 

P n + „ 
y Re 

oY y Re V 

. ay 

ax . 

J 


(2.49) 


CS 2 


respectively. In the above equations, the velocities have been 

nondimensionalized with respect to the characteristic velocities 

of the problems concerned, all lengths have been 

nondimensionalized with respect to the characteristic lenghts of 

the respective problems and the nondimens ional pressures have been 

2 

obtained with respect to Pi^char^ * 

We have applied the EXTRA-FLAG algorithm to five problems, 
namely, (i) Lid-driven square cavity (ii) Lid-driven oblique 
cavity (iii) Developing flow in a channel ( iv ) Investigation of 
confined wakes behind a square cylinder in a channel and (v) 
Investigation of confined wakes behind a circular cylinder in a 
channel. The characteristic velocity, length and the definition 
of Reynolds number are presented in Table 2.2 



Table 2.2 Characteristic Velocity, Length and Definitic 
Reynolds number for the Problems Under Consideration 


Problem Studied 

Characteristic 

Characteristic 


Velocity 

Length 

(i) Lid-driven square 

Lid velocity 

Height H 

cavity 

ULID 


(ii) Lid-driven oblique 

Lid velocity 

Oblique height 

cavity 

ULID 

H 

( iii )Developing flow in 

Average velocity 

Channel height 

a channel 

at the channel 

H 


inlet U 

av 


(iv) Investigation of 

Average velocity 

Cylinder height 

confined wakes 

at the channel 

B 

behind a square 

inlet U 

av 


cylinder in a 



channel 



fv) Investigation of 

Average velocity 

Cylinder 

confined wakes 

at the channel 

Diameter D 

behind a circular 

Inlet 


cylinder in a 



channel 
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The proposed EXTRA-FLAG algorithm has been embodied in a 
computer code to envisage the objectives enumerated above. The 
computations were done on the CONVEX C— 220 computer while pre-and 
post processing were done on a HP 9000 (M-340) Workstation. 

The results obtained by using the proposed algorithm have been 
discussed in the subsequent chapters. Code validation has been 
accomplished by comparing our results with available results in 
open literature. 
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Chapter-3 

3. Testing of the Algorithm on Model Problems 

The EXTRA— FLAG algorithm is tested on three flow 
configurations: driven cavity flow in square and oblique cavities 
and developing flow in a channel . In this section, the results 
obtained by the present algorithm are reported and further 
discussed how they compare with the results available in 
literature . 

3.1 Flow in a Lid-driven Square Cavity: 

The numerical domain consists of a two dimensional lid-driven 
square cavity with no-slip and impervious boundary conditions at 
the bottom and side walls except at the top , where u is a 
non-zero constant and equal to ULID. In the governing equations , 
the velocities have been nondimensionalized with respect to ULID. 
All lengths have been scaled with respect to cavity height H, and 
the Reynolds number is defined as Re = (tJLID.H)/^. A large number 
of numerical bench-mark solutions on this domain are available 

in open literature. 

The streamline patterns using a 41 x 41 grid and a 81 x 81 grid 
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for the aforesaid two dimensional cavity at Re = 400 are shown in 
Fig. 3. 1(a) and 3.1(b) respectively. The numerical solution of Ghia 
et al [110] compares favourably with the present solution 
calculated with 81 x 81 grid . As stated in the earlier chapterj 
a time-marching scheme was deployed and the flow became steady 
when the maximum deviation of each flow variable between 
consecutive time steps is below 5 x 10~^ . In Fig . 3 . 2 u velocity 
components at the vertical mid plane are plotted for 41 x 41 and 
81 X 81 grids and compared with the results obtained by Ghia et al 
[110] and cited by Peyret and Taylor [111] . It may be mentioned 
that Ghia et al employed a 129 x 129 grid via a multigrid 
technique . However, Fig. 3. 2 shows reasonably good agreement 
between Ghia et al ’ s result and the present computation with a 81 
X 81 grid . Even otherwise, the positions and values of the 
minimum u velocity in position and value of all the computations 
agree with each other on an overall assessment . The v velocity 
components at the horizontal mid plane of the cavity for the same 
Reynolds number (Re = 400) are shown in Fig. 3. 3, It is evident 
that the present computation with 81 x81 grid demonstrates 
extremely good agreement with that of Guj and Stella [112] for 
both the maximum and the minimum v velocity positions and values . 
However, the minor discrepancy with Ghia et al ’ s result ( with 
respect to maximum v velocity 7.5 percent, and with respect to 
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Extremes in Velocity for the Driven Cavity 
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3.2 Flow in Lid-driven Oblique Cavity 

Laminar flow in an oblique-cavity with moving lid at different 
Reynolds numbers is chosen as another important test case for 
which a non-orthogonal grid has been used. Fig. 3. 6(a) shows the 
grid and the boundary conditions. Fig. 3. 6(b) shows the velocity 
vectors in the cavity at a Reynolds number of 100 with a 21 x 21 
grid. The u-velocity distribution along the mid-plane for 
Reynolds numbers of 100, 400 and 1000 have been shown in Fig. 3. 7. 
The positions and values of minimum u-velocity for all the 
Reynolds numbers appear to be physically meaningful. Fig. 3. 8 
shows the velocity vectors in an oblique cavity at a Reynolds 
number of 1000. It is to be noted that our results for oblique 
cavity compare favourably with results of Peric [115]. Basically, 
these computations were done in order to test the ability of the 
algorithm in tackling non-orthogonal geometry. Due to sparse 
numerical or experimental data base on oblique cavity, rigorous 
comparison was not possible. However, this particular aspect of 
the algorithm has been tested on a complex geometry in a 
subsequent chapter. 

3 . 3 Developing Flow in a Channel 

Developing flow in a rectangular channel was calculated in an 


U= ULID, V = 0 




Fig. 3.7 Variation of U velocity along the mid plane 
in an oblique cavity. 



Fig. 3.8 Velocity vectors in an oblique cavity for 
Re = 1000 
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effort to find out the hydrodynamic entrance length in a two 
dimensional channel flow . In the governing equations ^ the 
velocities have been nondimensionalized with respect to the 
average incoming velocity at the channel inlet, all lengths 
have been nondimensionalized with respect to the channel height H 
and the Reynolds number is defined as Re = ( U^^.H ) / p . No slip 
and impervious boundary conditions for the axial and normal 
components of velocity are applied on the top and the bottom walls 
{u = 0, v = 0). At the inlet, the normal component of velocity is zero 
( v = 0 ) and a uniform axial velocity profile (u=U ) is 

a V 

deployed. At the exit of the channel, the second derivatives of the 
dependent variables in the flow direction are set equal to zero ( 
— - ~ — - = 0 ) in order to ensure smooth transition through the 

outflow boundary. The problem of laminar flow in the entrance 
region of ducts has been studied extensively by several 
investigators . We shall examine whether our results corroborate 
the results available in literature . 

Fig. 3. 9 shows the velocity profiles at the entrance region of a 
two dimensional rectangular channel for a Reynolds number of 50. 
In the immediate downstream of the inlet , the profile has a 
local minimum on the axis of symmetry and a pair of maxima located 
symmetrically on either side of the axis . In the farther 
downstream , the boundary layers grow in size ; consequently the 



so, C-^l-rL X Si gr-iczf) 


Fig. 3.9 Velocity profiles for devoloping flow in a two 
dimensional rectanguior channel. 


Fig. 3.10 Variation of CfRe in a two dimensional channel 
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velocity profiles tend to become fully developed and after a 
distance ( entrance length ) it remains unchanged . Convexity is 
induced in the velocity profile and finally it culminates as a 
parabolic profile . However , this trend of variation of velocity 
profile can be corroborated and expatiated by comparing the 
results of Abarbanel , Brandt et al [1161. 

Variation of skin friction ( C^Ee ) along the channel walls for 
different Reynolds numbers has been shown in Fig. 3.10. It may be 
mentioned that the distributions of skin friction ( C^Re ) along 
the channel walls are exactly same for the top and the bottom 
walls . As expected , the skin friction values ( C^Re ) for all 
the Reynolds numbers are very high near the inlet . Thereafter , 
they take a sharp plunge and asymptotically reach a constant value 
of 12 { which is the value corresponding to the fully developed 

Poiseuille flow ) . However , from our results it is quite evident 

that the entry length ( the length at which fully developed 
velocity profile culminates ) is a function of Reynolds number and 
this functional relationship may be summarized as 
(x / H ) « 0 , 055 Re 

It may also be mentioned that here we have considered 
attainment of the 99 % of velocity profile as the index of fully 
(jgveloped condition . If we consider that the entry length is made 
^p of inlet and filled regions i the estimation of entry length 
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will be somewhat more { see Mohanty and Das [117] ) . However > we 
are not estimating the value of entrance region to a very high 
degree of correctness in this analysis . Our purpose is to examine 
the overall validation of the proposed numerical algorithm . 
Within the framework of 99 % attainment of Poiseuille profile as 
definition of fully developed condition and corresponding entry 
length , our prediction of the same corroborates the results 
available in literature ( see Wang and Longwell [118] ) . 


Chapter 4 
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4. Study of Confined Wakes Behind a Square Cylinder in a Channel 


4.1 Introduction 

The oscillation of chimney stacks and other structures in 
transverse flows is caused by vortex shedding. An initially smooth 
and steady flow across a bluff body may bring about damaging 
deflection of the body in case the natural frequency of the 
obstacle is close to the shedding frequency of the vortices. If we 
concentrate on the wake region where Karman vortex street has been 
formed behind a bluff body, we shall observe that the wake zone 
undulates like a flag from side to side. This alternating 
deflections of the wake induce periodicity in the entire flow 
field. As a result the forces on the bluff body become periodic 
which culminate in vibration that can be detected from the 
oscillation of the bluff body. If this excitation __^frequency 
synchronizes with the natural frequency of the bluff body, the 
phenomenon of resonance is the obvious outcome. Hence the unsteady 
flows about the bluff bodies are of direct relevance to design of 
structures, road vehicles, heat exchangers and where ever the flow 
induced vibration is important. 

In this chapter, the numerical investigation of the structure 
of confined wakes behind a square cylinder is presented. Details 
of the vort ex~shedd ing phenomenon are simulated through numerical 
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flow visualization. We have used the well known MAC (Marker And 
Cell) algorithm of Harlow and Welch [1] to analyse the basic flow 
features of the problem. Subsequently the proposed EXTRA-FLAG 
algorithm has also been applied, in order to look into the similar 
aspects of the problem and compare the results with those obtained 
by MAC algorithm. 

It has been known from both experimental and numerical 
studies that at moderate Reynolds numbers vortex shedding behind 
the cylinder introduces periodicity in the flow field. The 
unsteady periodic wake can be characterized by the Strouhal number 
which varies with Reynolds number and blockage-ratio of the 
channel. The periodicity of the flow is, however, damped in the 
downstream of a long duct. This damping may be attributed to the 
influence of side walls on the flow structure. In the present 
study, such fundamental flow features as well as the 
non-dimensional releat ionships characterizing the vortex- shedding 
frequency, have been analysed in detail. 

4.2 Previous Work 

Vortex structure in the wake of a circular cylinder has been 
investigated both experimentally [119,120] and numerically 
[121,122]. These investigations have shown that vortex shedding 
behind a circular cylinder in an unbounded medium starts around a 
Reynolds number of 40 and periodicity is induced in the flow field 
[123]. Study of vortex shedding behind rectangular / square 
cylinders has also been the subject of Investigation for many 
researchers. A systematic study of eddies behind a rectangular 
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cylinder has been undertaken by Okajima [124] His experimental 
results show how Strouhal number varies with the aspect ratio of 
cylinders in the range of Reynolds number between 70 and 2 x 10^ A 
recent work of Okajima [125] presented the variation of lift and 
drag forces , base pressure and Strouhal number of rectangular 
cylinders with Reynolds number. His computations by finite 
difference method showed good agreement with experimental results 
and it was possible to detect a critical range of Reynolds number 
where the value of Strouhal number changes followed by a drastic 
change in flow pattern . Kelkar and Patankar [126] have computed 
steady two dimensional flow around a square cylinder at different 
Reynolds numbers and determined the onset of unsteadiness through 
a linear stability analysis of steady flow. Stability of the 
steady flow to small perturbations is examined by computing the 
evolution of these perturbations.Davis , Moore and Purtell [127] 
reported results of a numerical and experimental study of flow 
around a rectangular cylinder in a horizontal channel. Strouhal 
numbers obtained from their computations are in good agreement 
with their measurements. Baba and Miyata [128] have studied the 
vortical structures of the flow around a rectangular cylinder by 
numerical integration of Navier-Stokes equations and explained 
some features of nonlinear interaction between vortical motions of 
different scales. However, if the square cylinder is confined in a 
channel , irrespective of shedding of vortices in the near wake, a 
parabolic velocity profile will evolve again at the exit of a long 
channel [129]. Biswas, Laschefski, Mitra and Fiebig [129] studied 
structure of laminar wake and heat transfer in the presence of 
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thermal buoyancy m a horizontal channel with a built in square 
cylinder. Based on these investigations carried out on related 
topics, it can be said that the channel walls exert damping 
effects on periodic flow. But concrete inferences about the 
relationship between the wake-zone aerodynamics and the channel 
walls have not been drawn so far. The purpose of the present work 
is to perform numerical investigation of the influence of Reynolds 
number and the channel walls on the structure of wakes behind a 
square cylinder in a two dimensional duct. 


4.3 Statement of the problem 

The system of interest is a horizontal channel with an 
obstacle in the form of a square cylinder placed inside it (Fig. 
4,1). The dimensionless equations for continuity and momentum may 
be expressed in the following conservative form 


D 


au dv 

ax dY 


0 


(4.1) 


dv dU^ d(uv) _ dp , -i f 

dt ^ dX dY ~ ~ ax Re ' ^2 ^^2 ^ 


(4.2) 


du a(UV) dv^ _ dP , ^ t d^v ^ d^V ^ 
at dx ^ aY ~ ■ dY Re ^ 3^2 


(4.3) 


Boundary conditions of interest in this investigation are : 
At the top and bottom surfaces of the channel 
u = V = 0 , 

at the entrance to the channel 


t 
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u 

XJav 


1.5 



V = 0. 

At the exit of the channel, the continuity boundary condition is 
used by setting 

- 0 . 

ax^ 


This ensures smooth transition through the outlet. No-slip 
boundary conditions are used for the velocity components on the 
obstacle. The differential form, of governing equations presented 
above have been employed for solving the problem by the MAC 
procedure. While obtaining the predictions of EXTRA-FLAG for the 
sake of comparison, the integral equations (2.47-2.49) described 
in chapter 2 have been used. 


4.4 Method of Solution 

In order to solve the problem with EXTRA-FLAG algorithm the 
computational domain is divided into cartesian cells. The 
velocity nodes .are located at . the vertices and the pressures are 
located at the centre ot the cells. The continuity control 
volumes and the momentum control volumes are shown in Fig. 4. 2. 
The solution procedure through the EXTRA-FLAG algorithm has 
already been discussed in Chapter-2. This is not being elaborated 
here for the sake of brevity. For high Reynolds number situation. 


however, a minor modification has been introduced in the solution 
procedure for providing upwind bias. This modification in the 

algorithm is described below. 



H 



L 


in a horizontal channel with built-in obstacle 



Vessure • e Velocity 

I Grid layout in the comutational domain for 
EXTRA- FLAG algorithm 
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441 Implementation of Upwinding in EXTBA-FLAG 

In the study of flow around a square cylinder, it was 
observed that almost upto a cell Reynolds number of 6, no 
upwinding was needed for EXTRA-FLAG. However, for very high flow 
Reynolds numbers (with a cell Re of the order of 10 or more) , an 
upwind bias was necessary to obtain numerically stable results. 
In the present work, an elementary attempt has been made to 
introduce upwind bias in the convective terms (refer equations 
(2.25) and (2.26)). The scheme implemented is as follows. 
Considering the convective contribution to the x-momentum 

equation, it is redefined for upwinding purposes as 
[C^j] [u.] = ^"1^2 (u n^ + V ny)dl] [u.] (4.4a) 

where the interpolation functions are given extra bias in the 

flow direction. There are calculated by the expression 
9 9 

X W. u. = 5] (1 - a) N. u. + a u (4.4b) 

j^=l ^ j = l ^ ^ ^ 

where the subscript r corresponds to the node situated upstream of 
the concerned face. The weight oc lies between 0 and 1. 

Using the above approach, for instance, on the east face of 
the control volume (see, Fig.2.4), if the mass flow is outwards 
r=2 (central node) and if the mass flow is inwards r = 4 (east 
side mid-node). In a similar fashion, for all the faces of the 
control volume, upwind— weighted interpolation has been applied for 
both the x-momentum and 4-momentum equations. 

Admittedly, the above implementation is only a first step 
towards a scientific way of accounting for dominant convective 
effects. More elaborate study is required to develop proper 
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upwinding strategies for high Reynolds number flows. As mentioned 
earlier, we have also solved the problem with the MAC algorithm to 
validate the EXTRA-FLAG algorithm on one hand and to highlight 
some intricate flow physics of the problem on the other. In the 
subsequent section, we shall discuss in brief, the MAC algorithm. 


4.4.2 Solution Procedure by MAC Algorithm 

A modified version of Marker and Ceil (MAC) method [1»5] is 
used to obtain the numerical solution of equations (4.1) “ (4.3). 
The computational domain is divided into cartesian cells. 
Staggered grid arrangements are used, in which velocity components 
are defined at midpoints of the cell sides on which they are 
normal (Fig.4.3). The pressure is defined at the centre of the 

cell . 


To start a computational cycle, guess-fields for velocity and 
pressure are used. From these fields, corrected pressure and 
velocity fields are obtained by pressure and velocity iteration 
through continuity equation, as discussed in the following 
paragraph In detail. The corrected velocity and pressure fields 
are, in turn, used in the Navier-Stokes equations (4.2) and (4.3) 

to evaluate velocities tor the next time increment. 

In each cycle, solutions for velocities are obtained in 

two folds. First, the velocity components are advanced in time 


using the previous state of flow to calculate 
caused by convection, viscous stresses, and 
ever a time step of duration 5t. The explicit 
not necessarily lead to a velocity field 


the accelerations 
pressure gradients 
time increment may 
with a zero mass 
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divergence in each cell. In the subsequent second fold, 

adjustment of pressure and velocity is done by a iterative 
process in order to ensure mass conservation in each cell. This 
iterative correction of explicitly advanced velocity field through 
implicit continuity equation is equivalent to the solution of a 
Poisson equation for pressure [130]. The convective terms of the 
equations (4.2) and (4.3) are discretized by a weighted average of 
an upwind and central differencing scheme [21,131,132], Here, the 
formal first-order accuracy 0{AX) of the upwind differencing 
scheme retains something of the second order accuracy of the 

advection field [16]. Runchal and Wolfstein [133] used a similar 

kind of upwind scheme to compute driven cavity flows. Their 

results compare favorably with those due to second-order schemes. 

A concise description of the mathematical methodology for 
pressure-velocity iteration is described herein. The divergence of 
velocities D is calculated in each cell from equation (4.1). If 
the magnitude of D is greater than a prescribed small value £, the 
pressure in each cell is adjusted proportional to the negative of 
velocity divergency by 

5p? . = - u . (4.5a) 

J 1, j 

where w is defined as 

m =: m /{2at{-” + -^ )) (4.5b) 

® 5X 

The factor « is an over relaxation factor. After 5P has been 
added to the cell pressure, the velocities are consequently 



( 4 , 6 ) 
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+ (5t . 5P^^^)/5X 


yH + ^ 1 yH +^ 1 


- (6t . 6p" )/5x 

i I J 


' " ''?/ ' * . 6p-^.,/sy 


^n + 1 ^n + 1 

1,3-1 i*1'l 


( 6t . 6p^ . ) /6 y 

1 } J 


(4.7) 

(4.8) 

(4.9) 

(4.10) 


This scheme of solution is continued until a steady or periodic 
flow is obtained. The procedure automatically adjusts the 
pressure level during each time cycle and thus avoids the need to 
prescribe pressure boundary conditions during the solution of the 
Poisson equation for pressure [111]. 

In order to set the initial condition for tangential 
velocities, U. . at each cell is taken equal to unity, i.e.,^ — =1. 

Consequently, the transverse component of velocity V. . at each 

1 > J 

ceil is taken as zero. 

The conditions necessary to prevent numerical instabilities 

are determined from the Courant-Friedrichs-Lewy (CFL) condition 

and the restriction on grid-Fourier number . 

According to the Courant—Friedrichs—Lewy conditions , the 

distance the fluid travels in one time increment must be less than 

one space step. Therefore, 

r-— . r f A 1 1 N 

5t < min (4.11) 

When the viscous diffusion terms are more important, the 
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condition necessary to ensare stability is dictated by the 


restriction on the grid Fourier numbers, which results in 


pdt < 


1 (5x^ Sy^) 

^ (SX^ + Sy^) J 


<4.12) 


Final for each time increment is minimum of the 5t’s obtained 
from (4*11) and (4.12). 


4.5 Results and Discussion 

4,5.1 Selection of Numerical Grid and Computational Parameters: 

For these computations with MAC algorithm 200x34 and 396x66 

grids have been used . The computational results for 200 x 34 and 

398 X 66 grids show an average difference of 3% in the peak value 

of C„Re on the channel walls. However, the computation time with 
f B 

396 X 66 grids is nearly 9 times larger than with 200 x 34 grids. 
It was found that for all practical purposes 200 x 34 grids can 
produce grid independent results, although for some calculations 
396 X 66 grids were used. Uniform grids are deployed throughout 
the calculation domain. Computations have been carried out in a 
channel of length L / H = 6.125. The geometrical centre of the 
square cylinder Is located at a distance X„ = 2.125 from the 

inlet. The channel and cylinder axes are aligned < Fig. 4.1 ). The 
aspect ratio { A / B ) of the cylinder is 1 . The influence of 
different blockage ratios (B/H), namely, 0.125, 0.25, 0.3125 and 

0.375 on vortex-shedding was studied. Reynolds number was used as 

an input parameter. 
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4.5.2 Validation of Computed Results and Comparison Between the 
Predictions of MAC and EXTRA-FLAG 

The structure of the wake and its functional relationship 
with Reynolds number can be observed in (Fig. 4. 4) and (Fig. 4.5). 
In Fig. 4.4 a computed steady solution for B/H = 0.25 and Reynolds 
number 3 7 is shown. The recirculating wake is extended to 
nearly twice the obstacle width in the downstream, but it is 
steady (symmetrical) and vortex shedding has not started. At Re 

B 

- 85 the wake loses its original symmetry and the flow becomes 
periodic in the near wake ( Fig . 4 . 5 ) .The periodicity is suppressed 
in the downstream of the square cylinder by the channel walls and 
the flow at the exit of the channel tends towards steady 
parabolic. Davis, Moore and Purtell ( 127 ] observed periodicity in 
computations with the same geometrical configuration as ours, at 
Re^^ 100. Okajima [124] found periodicity in the wake behind a 
rectangular cylinder in an infinite medium ( blockage being 
negligibly small ) at Re = 70 . Shair et al [120] report a 

o 

critical Reynolds number of 130 for the wake to become periodic in 
their experiments with a circular cylinder in a channel. The 
blockage ratio was 0.33 and the Reynolds number was defined in 
terms of a maximum velocity U which would exist at the same 
location as that of the centre of the cylinder under flow 
conditions identical with those of the experiment but in absence 
of the cylinder. However, their transition Reynolds number 
describing the onset of periodicity, based on a definition similar 
to ours, would be 87 ( the average velocity is two-third of the 
maximum velocity ) . It can be said that our Reynolds number for 
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the onset of periodicity lies well within the range of critical 
Reynolds numbers obtained by other researchers- Finally, we may 
mention that the Figure 4.4 and 4.5 are not merely qualitative 
since the scales of the average channel velocity Uav are shown in 
the vector— plots and these can be used for quantitative analysis. 

Numerical calculations with higher Reynolds numbers (Re > 85) 
confirm shedding of vortices into the stream. As such, with 
increasing Reynolds number ( beyond the aforesaid transition 
Reynolds number ) , the von Karman vortex street is formed and 
alternate shedding of vortices into the stream becomes prominent. 
The vortex shedding and formation of the von Karman vortex street 
is better understood from the numerical flow visualization of Fig. 
4.6. For a Reynolds number of 162, separation is observed at the 
leading edge followed by rolling up of vortices behind the 
cylinder. The flow is seen to completely detach itself on the 
lower surface of the cylinder. A favourable comparison between the 
experimental and numerical flow-visualization has been reported by 
Okajima [124]. Our observation (Fig. 4.6) of numerical flow 
visualization follows a qualitative trend similar to Okajima s 
study [124] for a Reynolds number of 150. At Re^^ = 375 formation 
of von Karman vortex street and its serpentine bends become 
prominent (Fig. 4,7). Flow separates at the leading edge and does 
not reattach during a period of vortex shedding into the wake. 
The results discussed so far have been obtained with the help of 
the MAC algorithm. Similar results were obtained form our 
EXTRA-FLAG algorithm also. Next, we shall illustrate some of the 

the EXTRA-FLAG algorithm. Fig. 4. 8 shows the 


results due to 
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asymmetrxc wake and the onset of formation of the Von Karman 
Vortex street. Fig. 4. 9 shows a typical vortex street behind a 
cylinder for a Reynolds number of Re^ = 125. Fig. 4. 10 shows the 
velocity vectors in a channel where the entire wake zone and 
downstream undulates from side to side. However this undulation 
and the von Karman vortex street extends upto the end of the 
channel which is evident from the streamline plot (Fig. 4 .11). 

The proposed EXTRA-FLAG algorithm has not only shown the 
correct trend of the Wake— zone aerodynamics j it has also predicted 


quantitative values which are very close to the predictions by MAC 


algorithm. In the subsequent section, we shall discuss some more 


intricate physical aspects of partially confined wakes behind a 


square cylinder. The discrepancy between the quantitative values 


obtained from two different numerical methods, namely, MAC and 


EXTRA-FLAG are so insignificant that we are only highlighting the 


physical phenomena without mentioning the methods deployed for the 


results . 


4.5.3 Lift and Skin Friction Characteristics 


Fig. 4. 12, (a) and (b), show the time evolution of the lift 


coefficient for two different Reynolds numbers. In each case, 


periodic flow field is observed after a short transient flow. The 


vortex shedding frequency can be determined from the 
time— evolution plot of lift coefficient distribution. Time period 
T can be calculated computationally by observing the 
nondimensional time when the lift coefficient is just crossing the 
mean value. The difference between two such alternate time values 
( also shown in Fig. 4. 12) gives the time period T. Once the time 
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period IS known, the corresponding frequency ( f = 1 / T ) and the 
Strouhal number ( S = fB / Uav ) can be evaluated. The periodicity 
is characterized by the vortex shedding frequency, or so to speak, 
by the dimensionless Strouhal number. Variation of the Strouhal 
number over a large range of Reynolds number is discussed 
subsequently. 

Flow field in front of the cylinder seems to be nearly 
independent of the structure of the wake [134]. Influence of the 
location of the obstacle in the channel { ) is not of great 

significance with respect to wake structure. However, for a 
uniform entry-profile the effect of flow development will come 
into play. The effect of periodicity on shear stress distribution 
is shown in Fig. 4. 13 ( a and b ). It is evident that the local 
skin friction coefficient C { = C^Reg) on the walls at the channel 
entrance is equal to 12 ( the value of the fully developed laminar 
flow in a 2-0 channel ) and tends to this value far downstream of 
the cylinder. The deviation of skin friction coefficient from 12 
in front of the cylinder ( nearly at a distance of 2.5 B upstream 
of the cylinder ) shows the upstream influence of the obstacle. It 
is also seen that the distribution of the skin friction C ( 
C^Eeg ) is same on both the channel walls for a Reynolds number of 
Reg= 25 ( steady flow ). However, for unsteady periodic flow, the 
situation differs considerably. For a Reynolds number of 125, we 
can see some remarkably changed trend in the distribution of the 
skin friction coefficient on the channel walls at the rear of the 
cylinder. Where the upper plate has a very high value of shear 
stress ( at X « 3.38 ), a minimum value of shear stress is 



Axial position of the 
square cylinder 
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observed on the lower plate. However, the trend may be vice-versa 
at another location (say, X 3.73 ) . This oscillating structure 
IS damped gradually in the downstream of the cylinder. 

4.5.4 Frequency (Strouhal Number) Characteristics of Vortex 
Shedding 

Frequency (f) of vortex shedding was also measured by another 
technique. This technique borrows some idea from the measurements 
of Oka j ima [124] . We recorded the normal component of velocity at 
a position of 6 B behind the cylinder at y/H = 0.5 in the wake. 
For the steady flow, the flow smoothly divides and reunites around 
the cylinder. As a consequence, the normal component of velocity 
at the aforesaid location will be zero. However, for an unsteady 
periodic flow ( so to speak, for high Reynolds numbers ), the 
normal velocity component at the same location fluctuates. It is 
evident from Fig. 4.14(a) that the recorded signal of the 
fluctuating velocity in the wake results in a sinusoidal wave. The 
time period T can be calculated from such signal traces and the 
corresponding frequency f ( = 1/T ) and the Strouhal number S ( = 
fB/Uav ) can also be found out. Fig. 4.14(b) shows the oscillation 
of the lift coefficient for the same Reynolds number and 
geometrical conf iguration as those of Fig.4. 14(a) . Strouhal 

numbers obtained in both the cases completely agree with each 
other ( S = 0.164 ) . 

Effect of blockage ratio on the variation of Strouhal number 

with Re is shown in Fig. 4.15. With increasing blockage ratio the 
B 

value of Strouhal number increases. In all the cases, the Strouhal 
number undergoes a slight change with increasing Reynolds number. 





) I 234567891 

Non-dimensional time 
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Fig4.15 Effect of blockage ratio on variation of Strouhal 
number with Rgr 
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At a low Reynolds number, there is a steady reattachment 
behind the leading edge and the flow finally separates at the 
trailing edge which results in symmetrical vortices* At somewhat 
higher Reynolds number, after the separation at leading edge, the 
flow reattaches on either the upper or lower surface of the 
obstacle during alternate shedding of vortices into the stream. 
Further increase in Reynolds number, makes the separated flows 
detached from both the surfaces which eventually widens the wake 
accompanied by a relatively sharp change in Strouhal number. 
However, the range of Reynolds number at which the above-desired 
changes in the wake structure occur, depends upon the blockage 
ratio ( see Figs. 4.6, 4.7, 4.15 and Table 4.1). 

Table 4.1 also shows the Strouhal number calculations from 
two methods, namely, the signal traces of the fluctuating velocity 
components and the time-evolution of lift coefficients. The 
calculated values from both the methods agree with one another. 


TABLE 4.1: 

Strouhal number at different Reynolds number and Blockage Ratio 


Obser Re 
No. ® 


Strouhal Number ( S=fB / Uav ) 

from time-evolution from signal traces of 
of lift-coefficient fluctuating velocity.. 


1 

60 

0.125 

0.1521 

0.1520 

2 

70 

0.125 

0.1551 

0.1556 

3 

80 

0,125 

0.16 

0.159 

4 

100 

0.125 

0.1623 

0,162 

5 

125 

0.125 

0.164 

0.164 

6 

200 

0.125 

0. 1668 

0.1669 

7 

300 

0.125 

0.167 

0.167 
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Obser 

No. 

Re 

■n /u 

Strouhal Number 

( S=fB / Uav ) 


from time-evolution 
of lift-coefficient 

from signal traces of 
fluctuating velocity 

8 

400 

0.125 

0.1672 

0.1676 

9 

500 

0.125 

0.1672 

0.1676 

10 

600 

0.125 

0.1674 

0.1674 

11 

800 

0.125 

0.1651 

0.1654 

12 

87 

0.25 

0.236 

0.237 

13 

112 

0.25 

0.238 

0.238 

14 

150 

0.25 

0.2384 

0.2384 

15 

200 

0.25 

0.235 

0.235 

16 

250 

0.25 

0.212 

0.210 

17 

312 

0.25 

0.22 

0.22 

18 

375 

0.25 

0.212 

0.211 

19 

500 

0.25 

0.217 

0.216 

20 

625 

0.25 

0.213 

0.2128 

21 

87 

0.3125 

0.2845 

0.2853 

22 

112 

0.3125 

0.2875 

0.2887 

23 

150 

0.3125 

0.2913 

0.2897 

24 

200 

0.3125 

0.2914 

0.2886 

25 

250 

. 0.3125 

0.2885 

0.2859 

26 

312 

0.3125 

0.2842 

0.2778 

27 

87 

0.375 

0.3317 

0.3316 

28 

112 

0.375 

0.3452 ■ 

0.3454 

29 

150 

0.375 

0.3488 

0.3489 

30 

200 

0.375 

0.3506 

0.3511 

31 

250 

0.375 

0.3566 

0.3423 

32 

312 

0.375 

0.3337 

0.3203 


Fig. 4. 16 compares the Strouhal number distribution with the 
experimental results of Okajima [124]. Strouhal numbers show 
continuous but slight change with Reynolds number . The present 
solution yields a somewhat higher Strouhal number for the entire 


© Present computation 
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Fig 4.16 Variation of Strouhal number with Reynolds number 
for a square cylinder 
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range of Reynolds numbers. This discrepancy between the 
experimental and calculated results can be explained by noting 
that the experiments were done for a negligibly small blockage 
ratio. Although our predictions were done for a very small 
blockage ratio ( B/H - 0.125 ), one may conjecture that the 

influence of side walls could not be completely ignored. The 
finite blockage might have brought about some change in Strouhal 
number . 

The predicted fluctuating velocity signals were monitored at 
various other locations in the duct. The fre<g.uencies of 
oscillation for a particular situation ( for a fixed Reynolds 
number and blockage ratio ) at different spatial locations were 
found to agree with one another. In the far downstream the 
amplitude becomes so low that the oscillation becomes 
insignificant. Possibly because of this reason, in a long duct, 
von Karman vortex street is gradually damped out and a steady 
parabolic profile reappears near the exit. 
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Chapter- 5 

5» Study of Confined Wakes Behind a Circuiar Cylinder in a 

Channel 

5,1 Introduction 

The vortex structure in the wake of a circular cylinder in an 
infinite medium has been the subject of numerous theoretical 
investigations. These studies include semi-analytic solutions of 
Navier-Stokes equations for flow past a cylinder due to Patel 
1 135, 136], numerical studies of Thoiaan and Szevrczyk [121], Lin et 
al [137], Fornberg [138] and Tuann and Olson [139]. Experimental 
investigations have also been undertaken to identify the structure 
of wake behind circular cylinder for different ranges of Reynolds 
number by Acrivos et al [119], Shair et al [120] and West and 
Apelt [140], Coutanceau and Bouard [141,142] used a visualization 
method to obtain the main, features of the flow field behind a 
circular cylinder in a tank. The change in the geometrical 
parameters describing the eddies with Reynolds number as well as 
the ratio between the diameters of the cylinder and the tank are 
presented. Their study also predicts the range of Reynolds 
numbers for which a closed wake exists and adheres stably to the 
cylinder. Braza et al [143] have computed the wake - zone — 
aerodynamics for incompressible flow past a circular cylinder 
using a second order accurate velocity-pressure formulation. They 
have discussed the interaction between the pressure and velocity 
fields over a range of Reynolds numbers, 100 ^ Re ^ 1000. The 
method used is a predictor-corrector pressure scheme, similar to 



that proposed by Amsden and Harlow [2] in SMAC. The vortex 
shedding was initiated by applying a physical perturbation imposed 
numerically for a short time. However, it has been recognized 
that there is very little information on the wake-structure when 
the cylinder is partially enclosed by side walls. Taneda et al 
[144] presented some observations of the flow past a circular 
cylinder towed through water near a plane wall. For a very small 
gap between the wall and the cylinder, they observed a single row 
of vortices are shed at a Reynolds Number of 170. Kiehm, Mitra 
and Fiebig [134] have solved full two and three dimensional 
Navier-Stokes equations behind a circular cylinder in a channel at 
different Reynolds numbers and blockage ratios. They have found 
the development of helical vortex tubes in the wake of a circular 
cylinder for three dimensional flows. 

However, the primary purpose of our study is to see how our 
EXTRA-FLAG algorithm performs in a multiply connected complex 
geometry involving a rectangular channel with a circular cylinder 
kept in the middle. We also aim at determining the influence of 
channel walls on the wake flow behind the cylinder. 

5.2.1 Flow Model 

From the knowledge of flow past a cylinder in an infinite 
medium, it can be said that the existence of a wake bubble 
(Fig. 5.1) behind the cylinder is usually observed for Re>5. The 
length of the wake bubble grows with increasing Reynolds 
number- In an infinite medium, beyond a Reynolds number of 40, the 
wake loses it symmetry and vortex shedding starts. Although 
information on partially confined wake is limited, Shair et al 
[120] and Kiehm et al [134] have shown that the critical Reynolds 
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number at which periodicitv . 

y arts, grows with increasing blockage 

r^tio (which is defined by the ■rat-T.-v j. 

oy nne ratio of diameter of the cylinder 

to height of the channel). 

Arising out of the above mentioned observations and the 
existing knowledge from open literature we recognise that the flow 
field around a circular cylinder in a channel is characteised by 
the following parameters 

the Reynolds number 

the blockage ratio (D/H), and 

the velocity profile at the channel inlet. 

Consequently , a computational domain which consists of two 
solid walls with no-slip conditions and a circular cylinder 
between these walls, is selected as shown in Fig. 5. 2. No-slip and 
impervious boundary conditions are applied on the whole surface of 
the circular cylinder. At the inlet of the channel, parabolic 
profile for the axial velocity has been deployed and the 
transverse component of the velocity has been taken as zero. At 
the exit of the channel, second derivatives of the velocity 
components in the flow direction have been set equal to zero. 
This ensures smooth transition of flow through the outflow 
boundary. 

5.2.2 Basic Equations 

The basic equations to be solved for the calculation of the 
laminar incompressible flow are the mass-conservation and the 
time-dependent Navier-Stokes equations with primitive variables. 
In nondimens ional form these are expressed as equations (2.47), 
(2.48) and (2.49) [see chapter-2]. The Reynolds number is defined 
by Re = P- where U stands for the average velocity at the 

JJ ^ filV 
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channel inlet and D denotes the cylinder diameter, 

5,2,3 Method of Solution 

EXTRA-FLAG algorithm is employed to obtain numerical 
solutions of the flow equations. For the implementation of 
numerical solution, the computational domain (Fig. 5. 2) is divided 
into a number of cells (see Fig.5.3(a)). Fig.5.3Cb) shows an 
enlarged portion of the mesh highlighting the zone surrounding 
the cylinder. A (99x33 ) grid is used for computation. The velocity 
components are collocated and they are defined at the vertices of 
the cells while the pressure has been defined at the centre of the 
cell. The details of solution procedure have been outlined in 
chapter-2. Upwinding was used in the evaluation of convection 
terms. However, our way of handling convective kinematics has 
already been mentioned in Chapter-4, We are not mentioning 
anything more on the algorithm for brevity. 


5.3 Results and discussion 


Computations were performed at Reynolds numbers of 87, 112 

and 625 for the blockage ratio of 0.25. It may be mentioned that 
for none of the above mentioned cases steady state solution could 
be obtained. However, the solutions were identified as unsteady 
periodic. Periodloty was characterized by the vortex shedding 
frequency which was obtained from the signal tracing of the 
transverse velocity components at a location 6D behind the 
cylinder on the centre-line of the channel. As pointed out 
chapter-4, this technique borrows the idea from the experimental 
frequency-determination technique of Okajima [124]. However, for 
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the first two circumstances (Re^^ - 87 and 122), the Strouhal 
numbers (non-dimensional frequency) were very close to each other 
with a magnitude of 0.2341 while for = 625 the Strouhal number 
was 0.2409. 

Fig. 5. 4 shows the velocity vectors in a channel with built-in 
circular cylinder for a Reynolds number of 122. The flow has 
become periodic in the near wake. The periodicity is suppressed 
in the downstream of the cylinder by the channel walls and the 
flow at the exit of the channel again assumes an almost parabolic 
profile . 

Fig. 5. 5 shows the streamlines crossing the cylinder in a duct 
at a Reynolds number 112. Alternate shedding of vortices and 
initiation of von Karman vortex street is quite evident from the 
figure . 

Fig. 5. 6 shows the streamlines in a channel with built-in 
circular cylinder for a Reynolds number 625. Formation of von 
Karman vortex street is clearly seen in this figure. With an 
increase in Reynolds number, the periodicity of the wake becomes 
more prominent. It can also be appreciated that despite the 
damping influence of side walls, the undulation of the wake is not 
completely suppressed in the downstream. 

However, the primary purpose of taking up this problem was to 
test the capability of EXTRA— FLAG algorithm. After analyzing the 
presented results, possibly it can be said that in the range of 
Reynolds number studied in this problem, the algorithm works 
satisfactorily. The study of wake structure for flow past a 
circular cylinder is indeed a difficult problem. As far as flow 
physics 13 concerned, we do not intend to draw any conclusion fro 
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Velocity vectors in a channel with built in 
circular cylinder, Rep^ 11^ 
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our brief study. The formidable task of the detailed simulation 

of this complex problem may be identified as the scope for future 

extension. 
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APPEND X A 

Table for CPU time requ red on COMVEX c 220 mach ne for 
“'rent flow problems dea t w th hi chapter 3 4 and 5 


PROBLEM 

GRID 

TIME MONITORED FOB 

Square Cavity Problem 
for Re = 10 0 

21 X 21 

steady state solr 

square cavity problem 

21 X 21 

steady state so In 

for Re = 1000 



Oblique cavity problem 

21 X 21 

steady state so In 

for Re = 1000 



Developing flow in a 

21 X 11 

steady state so In 

channel for Re = 40 



♦Flow over square cylin- 
der placed in a channel 
solved by EXTRA FLAG 
for Re , = 45 0 

97 X 33, 

1 

1 

on e t i mes t e p ca 1 - 
culation after 
s t ab I 1 i za t 1 0 n of 
periodicity 

•Flow over square cylin- 
der placed in a channel 
solved by M_A C method 
for Re != 450 

200 x34 

one times top cal- 
culation after 
s t ab i I i za t i on of 
period ici ty 

•Flow over circular cy 1 in 
der placed in a channel 
solved by EXTRA FLAG 
for Re = 450 

97 X 33 

one t t mes t ep ca 1 - 
cu lat i on af ter 
stabilization of 
period ici ty 


The current implementation of the EXTRA FLAG may lack in 
tization in terms of efficient use of skilled programming. The 
time and storage requirements are parameters which are very 
tive to the particuiar implementation of the algorithm, 
over, the transient problems fail in the category of impulsive 
which introduces spatial singularities in continuum or 
sntrated imbalance of mass. Hence the first few time steps may 
up to several hundreds of ’continuity / pressure-velocity 
Ltions However, after a few time steps . the mass errors 

diffused in space and the number of continuity iterations goes 
1 to a few tens or even less. 
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